```{r}
# --- Load required packages ---
library(readxl)
library(dplyr)
library(lmtest)
library(sandwich)
library(stargazer)

# 1) Load the data
df <- read_excel("Data Model 3A and 3B.xlsx",
                 sheet = "Planilha1")

# 2) Rename key columns
df <- df %>%
  rename(
    Employment = Yit_Emprego,
    TAX        = `NSit_Nao Simples`,
    Post       = POSLEIt_Post
  )

# 3) Create PAY dummy from 'Situacao'
df <- df %>%
  mutate(PAY = as.integer(Situacao == "Reonerado"))

# 4) Build interaction variables
df <- df %>%
  mutate(
    TAXxPAY      = TAX * PAY,
    TAXxPost     = TAX * Post,
    PAYxPost     = PAY * Post,
    TAXxPAYxPost = TAX * PAY * Post
  )

# 5) Define RHS variable vectors
base_vars <- c("TAX", "PAY", "Post",
               "TAXxPAY", "TAXxPost", "PAYxPost", "TAXxPAYxPost")

div_vars <- c("DummyDivCNAE_02","DummyDivCNAE_03","DummyDivCNAE_08","DummyDivCNAE_10",
              "DummyDivCNAE_13","DummyDivCNAE_14","DummyDivCNAE_15","DummyDivCNAE_17",
              "DummyDivCNAE_18","DummyDivCNAE_20","DummyDivCNAE_21","DummyDivCNAE_22",
              "DummyDivCNAE_23","DummyDivCNAE_24","DummyDivCNAE_25","DummyDivCNAE_26",
              "DummyDivCNAE_27","DummyDivCNAE_28","DummyDivCNAE_29","DummyDivCNAE_30",
              "DummyDivCNAE_31","DummyDivCNAE_32","DummyDivCNAE_33","DummyDivCNAE_41",
              "DummyDivCNAE_42","DummyDivCNAE_43","DummyDivCNAE_47","DummyDivCNAE_49",
              "DummyDivCNAE_50","DummyDivCNAE_51","DummyDivCNAE_52","DummyDivCNAE_55",
              "DummyDivCNAE_58","DummyDivCNAE_60","DummyDivCNAE_62","DummyDivCNAE_63",
              "DummyDivCNAE_82")
year_vars <- c("DummyAno_2014","DummyAno_2015","DummyAno_2016","DummyAno_2017",
               "DummyAno_2019","DummyAno_2020","DummyAno_2021")

all_vars <- c(base_vars, div_vars, year_vars)

# 6) Build the formula and fit Model 3A (Employment as DV)
fmla <- as.formula(paste("Employment ~ -1 +", paste(all_vars, collapse = " + ")))
mod3A_sel <- lm(fmla, data = df)

# 7) Print coefficient summary with robust SEs (HC1)
coeftest(mod3A_sel, vcov = vcovHC(mod3A_sel, type = "HC1"))

# 8) Stargazer summary table: HTML & LaTeX (with all FE)
stargazer(mod3A_sel,
          title            = "Table: Triple-Difference Effect on Employment (3A partial)",
          type             = "html",
          out              = "3A_partial_stargazer.html",
          digits           = 2,
          covariate.labels = c("TAX", "PAY", "Post",
                               "TAX×PAY", "TAX×Post", "PAY×Post", "TAX×PAY×Post"),
          omit.stat        = c("f", "ser"),
          header           = FALSE
)

stargazer(mod3A_sel,
          title            = "Table: Triple-Difference Effect on Employment (3A partial)",
          type             = "latex",
          out              = "3A_partial_stargazer.tex",
          digits           = 2,
          covariate.labels = c("TAX", "PAY", "Post",
                               "TAX×PAY", "TAX×Post", "PAY×Post", "TAX×PAY×Post"),
          omit.stat        = c("f", "ser"),
          header           = FALSE
)

message("Tables saved as 3A_partial_stargazer.html and .tex")

# Stargazer summary, omitting CNAE and year FEs, with explanatory note
stargazer(mod3A_sel,
          title             = "Model 3A: Triple Differences Estimator – Employment",
          type              = "html",
          out               = "3A_partial_stargazer.html",
          digits            = 2,
          covariate.labels  = c("TAX", "PAY", "Post",
                                "TAX×PAY", "TAX×Post", "PAY×Post", "TAX×PAY×Post"),
          omit              = c("DummyDivCNAE", "DummyAno"),
          omit.stat         = c("f", "ser"),
          add.lines         = list(
            c("Fixed effects", "CNAE divisions and years")
          ),
          header            = FALSE
)

stargazer(mod3A_sel,
          title             = "Model 3A: Triple Differences Estimator – Employment",
          type              = "latex",
          out               = "3A_partial_stargazer.tex",
          digits            = 2,
          covariate.labels  = c("TAX", "PAY", "Post",
                                "TAX×PAY", "TAX×Post", "PAY×Post", "TAX×PAY×Post"),
          omit              = c("DummyDivCNAE", "DummyAno"),
          omit.stat         = c("f", "ser"),
          add.lines         = list(
            c("Fixed effects", "CNAE divisions and years")
          ),
          header            = FALSE
)

message("Tables saved as 3A_partial_stargazer.html and .tex")

# 9) Subset to Employment > 0 for log outcome
df_pos <- df %>%
  filter(Employment > 0)

# 10) Create the log-dependent variable
df_pos <- df_pos %>%
  mutate(logEmp = log(Employment))

# 11) Fit Model 3B (logEmp as DV, same covariates, -1 for no intercept)
fmla3B <- as.formula(paste("logEmp ~ -1 +", paste(all_vars, collapse = " + ")))
mod3B_sel <- lm(fmla3B, data = df_pos)

# 12) Print summary with robust SEs (HC1)
coeftest(mod3B_sel, vcov = vcovHC(mod3B_sel, type = "HC1"))

# 13) Stargazer summary for log(Employment), omitting FE
stargazer(mod3B_sel,
          title             = "Model 3B: Triple Differences Estimator – Log(Employment)",
          type              = "html",
          out               = "3B_partial_stargazer.html",
          digits            = 2,
          covariate.labels  = c("TAX", "PAY", "Post",
                                "TAX×PAY", "TAX×Post", "PAY×Post", "TAX×PAY×Post"),
          omit              = c("DummyDivCNAE", "DummyAno"),
          omit.stat         = c("f", "ser"),
          add.lines         = list(
            c("Fixed effects", "CNAE divisions and years")
          ),
          header            = FALSE
)

stargazer(mod3B_sel,
          title             = "Model 3B: Triple Differences Estimator – Log(Employment)",
          type              = "latex",
          out               = "3B_partial_stargazer.tex",
          digits            = 2,
          covariate.labels  = c("TAX", "PAY", "Post",
                                "TAX×PAY", "TAX×Post", "PAY×Post", "TAX×PAY×Post"),
          omit              = c("DummyDivCNAE", "DummyAno"),
          omit.stat         = c("f", "ser"),
          add.lines         = list(
            c("Fixed effects", "CNAE divisions and years")
          ),
          header            = FALSE
)

message("Tables saved as 3B_partial_stargazer.html and .tex")
```

